rm(list = ls(all = TRUE))
setwd("C:\\Users\\lenovo\\Desktop\\R")
#10.1 相似性的度量
#10.1.1 样品间的距离
#欧氏距离、绝对值距离、切比雪夫距离、明氏距离、兰氏距离
#默认并使用最多的是欧氏距离
#method = euclidean
set.seed(1234)
dat <- matrix(rnorm(10), nrow = 5)
dat
dist(dat)
#距离矩阵为对称矩阵，dist默认输出只显示距离矩阵的下三角
dist(dat, diag = TRUE, upper = TRUE)

#变量的量纲不同、测量值变异程度相差悬殊时，先对数据进行标准化
set.seed(123)
dat <- matrix(c(rnorm(5), rnorm(5,mean = 100,sd = 10)),nrow = 5)
dat
colMeans(dat)
#第一列均值为0，第二列均值为100，第二列的权重比第一列大
dist(dat)
#使用函数scale将两列标准化后再计算距离矩阵
dat.scale <- scale(dat)
dist(dat.scale)

#有时，变量间有多重相关性，其造成的信息重叠，
#会片面强调某些变量的重要性
#改进方法：马氏距离
#利用函数mahalanobis计算样本点到某个中心点的马氏距离
#自主编写：
dist.ma <- function(data){
  #将数据化为矩阵
  x <- as.matrix(na.omit(data))
  #cov函数计算的是列与列的协方差
  V<- cov(x)
  #对于正定矩阵A，可对其进行Choleskey分解
  #Cholesky 分解，t为行列转换
  L <- t(chol(V))
  #用于解特殊线性方程组，系数矩阵为下三角
  stdX <- t(forwardsolve(L,t(x)))
  #算距离，存为矩阵
  d <- dist(stdX)
  #给输出的矩阵把标签加上
  attr(d, "Labels") <- row.names(data)
  #输出矩阵
  d
}

#可用于计算样本点之间的马氏距离
#其中的协方差矩阵用样本协方差矩阵估计
dist.ma(dat)

#10.1.2 变量间的相似系数
#如果要对变量进行聚类，需要度量变量之间的相似程度
#最常用的度量变量之间的相似程度的指标有相关系数和夹角余弦
#cor
set.seed(1234)
x <- rnorm(20)
y <- rnorm(20)
cor(x,y)
#夹角余弦计算
sum(x*y)/sqrt(sum(x^2)*sum(y^2))

#若将原始数据标准化，此时的样本相关系数与夹角余弦等价
x.scale <- scale(x)
y.scale <- scale(y)
sum(x.scale * y.scale)/sqrt(sum(x.scale^2)*sum(y.scale^2))

#10.2 层次聚类法
#系统聚类法，先分成几类，再将相似的两类合在一起，直到合并为一类
#合并类的过程中，需要度量类之间的相似性

#10.2.1 类之间相似系数的定义

#使用函数helust
#参数：method=
#single 最短距离法
#complete 最长距离法（默认）
#median 中间距离法
#average 类平均法
#centroid 重心法
#mcquitty 相似法
#ward.D2 离差平方和法


#10.2.2 Q型聚类
#加载数据
#install.packages("flexclust")
library(flexclust)
data(nutrient)
head(nutrient)

#将其行名（样品标签）换成小写字母
row.names(nutrient) <- tolower(row.names(nutrient))
head(nutrient)

#Q型聚类
#先查看变量的类型和分布情况
str(nutrient)
summary(nutrient)
#虽然变量皆为数值型，但是各变量量纲差距较大
#标准化
nutrient.scaled <- scale(nutrient)
#欧氏距离计算距离矩阵
d.eu <- dist(nutrient.scaled, method = "euclidean")
#层次聚类
hcl <- hclust(d.eu, method = "average")
#树状图绘制
#plot为泛用函数，实际上调用的是plot.hclust
#hang默认为正，换为-1，使得聚类图中样本的标签处于同一水平线上
plot(hcl, hang = -1)
#如果确定了类别的个数，
#可以使用函数rect.hclust在树状图的分支周围绘制矩形
#以突出相应的簇
rect.hclust(hcl, k = 5)

#一般来说，分析者可以根据实际问题确定需要的聚类个数
#也可以借助Nbclust包
#install.packages("NbClust")
library(NbClust)
NbClust(nutrient.scaled, distance = "euclidean", method = "average")
#因为函数NbClust将绘图画布设置为一行两列，用par将之回复为一行一列
par(mfrow = c(1,1))

#输出D index法的结果
#左图中下降最快的点（肘点）即为最优聚类个数,
#从4类到5类D index值下降得最快，之后很慢，因此选5类
#右图中D index值的二阶差分上升最快的点就是对应的最优聚类个数
#据此，聚类个数仍选为5

#输出的第二部分为通过改变聚类个数、距离度量方法和类之间距离度量方法
#得到所有组合下的聚类结果
#并据此向用户提出最优聚类个数的建议
#结果表明，建议聚类个数为2、3、5、15的都有最多的4个组合

#假定根据上述结果，选择5类
#显示每种食品所属的类别
cutree(hcl, k = 5)
#将上面的结果排序后放入一个数据框中
data.frame(group = sort(cutree(hcl, k = 5)))

#换一种聚类方法
#用马氏距离计算样品间的距离，并用离差平方和法进行层次聚类
d.ma <- dist.ma(nutrient)
hc2 <- hclust(d.ma, method = "ward.D2")
plot(hc2, hang = -1)
rect.hclust(hc2, k = 5)
data.frame(group = sort(cutree(hc2, k = 5)))
#不同的聚类方法结果会有差异，建议多方法尝试取较稳定的那种


#层次聚类法还可以用于混合类型的数据
#install.packages("cluster")
library(cluster)
data("flower")
str(flower)
summary(flower)
#从取值范围来看，flower数据的第7、8组数据取值范围不同，需要标准化
flower[,c("V7","V8")] <- scale(flower[,c("V7","V8")])
summary(flower)
#计算距离矩阵。因为数据框里有因子型变量，因此不用dist
#用cluster包里的daisy
#该包支持混合型变量的距离计算，
#非数值型变量采用高氏（Gower）距离度量变量之间的距离
dmat <- as.matrix(daisy(flower))
#显示前6行，前6列
dmat[1:6,1:6]
#用函数agnes进行层次聚类
#diss = TRUE 表明输入的是距离矩阵而非原始数据框，
#method取默认值average，表示类平均法
flower.cluster <- agnes(dmat, diss = TRUE)
pltree(flower.cluster, hang = -1, main = "")
#画出的图，横坐标表示样品的编号，纵坐标表示距离
rect.hclust(flower.cluster, k = 3)
#显示各个样品所属的类别
cutree(flower.cluster, k = 3)

#10.2.3 R型聚类

#R型聚类与Q型聚类本质上并无区别，
#只是将聚类的对象由样品换成了变量
#因此相似性度量从样本间的距离换成了变量
#即观察指标之间的相关系数

#输入数据
R <- matrix(c(1,0.852,0.671,0.099,0.234,0.376,
              0.852,1,0.636,0.055,0.174,0.321,
              0.671,0.636,1,0.153,0.233,0.252,
              0.099,0.055,0.153,1,0.732,0.637,
              0.234,0.174,0.233,0.732,1,0.676,
              0.376,0.321,0.252,0.627,0.676,1),
            nrow = 6,
            dimnames = list(c("身高","下肢长","手臂长","腰围","胸围","臀围")))
#两变量之间的相关系数的绝对值（不是相关系数）越大，则两变量越相似
#或称之为，距离越近
#1-|R|才是距离
#可以用abs函数将负的想关系数转换为正数
d <- as.dist(1 - R)
hc <- hclust(d)
plot(hc, hang = -1)
rect.hclust(hc, k = 2)
#大体可分为两大类，一类为反映人身高的变量，
#另一类为反映人胖瘦（或围度）的变量

#10.3 k均值聚类
#层次聚类法每一步都需要计算类之间的距离，样本过大时运算量过大
#动态聚类法：k均值聚类(k-means clustering)
#首先指定一个参数k，随机选取k个样品作为中心点，
#其余样品根据它们中心点的距离逐一归为其中一类
#同时计算新的中心点，即均值向量
#重复上一步，直到所有样品全部归类
#重复上述步骤，直到前后两次调整的结果一致

#设定中心点的个数为5
set.seed(1234)
km <- kmeans(nutrient.scaled, centers = 5)
km
#括号中显示组间平方和占总平方和的78.3%，类似于方差分析的思想，
#此百分比越大表明组内差距越小、组间差距越大，即聚类效果越好

#对聚类结果进行排序以查看分类情况
sort(km$cluster)

#k-means算法每次都是随机选取k个样品作为中心点，再计算距离，然后重新分类
#每次的运行结果都有差异
#对类别及k的数值选择不同，分类的结果也不完全相同
#选择一个正确的聚类数目对于划分样品比较重要
#k值可以根据专业知识确定，也可以先用层次聚类法确定个数


#使用fpc包里的函数plotcluster把聚类结果投影到两个区分度最高的维度
#以验证聚类的结果
#install.packages("fpc")
library(fpc)
plotcluster(nutrient.scaled, km$cluster)
#从图可知，聚类效果较好，如果效果不好，需要调整参数centers的取值
#或者换用别的聚类方法

#k-means聚类算法简单、聚类速度快，对于球形数据聚类的效果较好
#但k-means算法无法识别非球形的簇，可采用基于密度的聚类方法
#fpc包里的函数dbscan可以实现密度聚类

#cluster包提供了大量的聚类方法和可视化技术
#agnes可实现凝聚的层次聚类算法（如agglomerative）
#diana可实现分类的层次聚类算法（如divisive analysis）
#fanny可实现模糊聚类分析（fuzzy cluster analysis）
